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Abstract 

In previous work (Olshausen & Field 1996), an algorithm was described for learning linear sparse codes 
which, when trained on natural images, produces a set of basis functions that are spatially localized, 
oriented, and bandpass (i.e., wavelet-like). This note shows how the algorithm may be interpreted within a 
maximum-likelihood framework. Several useful insights emerge from this connection: it makes explicit the 
relation to statistical independence (i.e., factorial coding), it shows a formal relationship to the algorithm 
of Bell and Sejnowski (1995), and it suggests how to adapt parameters that were previously fixed. 



Copyright © Massachusetts Institute of Technology, 1996 

This report describes research done within the Center for Biological and Computational Learning in the Department of Brain 
and Cognitive Sciences at the Massachusetts Institute of Technology. This research is sponsored by an Individual National 
Research Service Award to B.A.O. (NIMH F32-MH11062) and by a grant from the National Science Foundation under contract 
ASC-9217041 (this award includes funds from ARPA provided under the HPCC program) to CBCL. 



1 Introduction 
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such that 



There has been much interest in recent years in unsuper- 
vised learning algorithms for finding efficient representa- 
tions of data. Among these are algorithms for sparse or 
minimum entropy coding (Foldiak 1990; Zemel 1993; 01- 
shausen & Field 1996; Harpur & Prager 1996), indepen- 
dent component analysis (Comon 1994; Bell & Sejnowski 
1995; Amari et al. 1995; Pearlmutter & Parra 1996), and 
hierarchical generative modeling (Dayan 1995; Hinton 
et al. 1995). One finds common threads among many 
of these techniques, and this note is an attempt to tie 
some of them together. In particular, I will focus on the 
sparse coding algorithm of Olshausen and Field (1996) 
and its relation to maximum-likelihood techniques. As 
we shall see, forming this link enables one to see a for- 
mal relationship to the independent component analysis 
algorithm of Bell and Sejnowski (1995), which although 
not originally described in terms of maximum-likelihood 
may be understood in this light. I shall also show how 
the algorithm may be cast in terms of mean-field theory 
techniques in order to obtain a lower bound on the log- 
likelihood, which shares some similarity to the use of a 
"recognition distribution" in the Helmholtz machine of 
Dayan et al. What emerges from this process is a bet- 
ter understanding of the algorithm and how it may be 
improved. 

2 Learning linear sparse codes 

In the sparse coding learning algorithm of Olshausen and 
Field (1996), a set of basis functions, <j>i{x), is sought 
such that when an image, I(x), is linearly decomposed 
via these basis functions, 



P\x) = y^a^ 8 -(ic) 



(1) 



the resulting coefficient values, a;, are rarely active (non- 
zero). In other words, the probability distribution over 
the a j should be unimodal and peaked at zero with heavy 
tails (positive kurtosis). This is accomplished by con- 
structing an energy function of the form 



E(I,a\</>) 
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(2) 
and then minimizing it with respect to the a; and <f>i. 
The first term in Equation 2 ensures that information is 
preserved (i.e., that the <f>i span the input space), while 
the second term incurs a penalty on activity so as to 
encourage sparseness. The intuition behind the choice 
of S is that it should favor among activity states with 
equal variance (|a| 2 ) those with the fewest number of 
non-zero (or not-close-to-zero) components. The choices 
experimented with include |a 8 |, log(l + a?), and — e a > . 

Gradient descent on E is performed in two phases, one 
nested inside the other: For each image presentation, E 
is minimized with respect to the a;; the <f>i then evolve 
by gradient descent on E averaged over many image pre- 
sentations. Stated more formally, we seek a set of basis 



arg min (min E(I,a\<f>) 

<p \ a 



(3) 



where (•) denotes an ensemble average over the images. 
Note that in this expression and in the rest that follow, 
I refers to the vector with components I(xj), a refers 
to the vector with components a;, and <f> refers to the 
matrix with components 4h(xj). 

The intuition behind the algorithm is that on each im- 
age presentation, the gradient of S "sparsifies" the activ- 
ity on the a j by differentially reducing the value of low- 
activity coefficients more than high-activity coefficients. 
This weeds out the low-activity units. The <f>i then learn 
on the error induced by this sparsification process. The 
result is a set of <f>i that can tolerate sparsification with 
minimum mean-square reconstruction error. A virtu- 
ally identical algorithm was developed independently by 
Harpur and Prager (1996). 

3 Maximum-likelihood framework 

While the energy function framework provides a useful, 
intuitive way of formulating the sparse coding problem, 
a probabilistic approach could provide a more general 
framework. Harpur and Prager (1996) point out that 
the first term on the right-hand side of Equation 2 may 
be interpreted as the negative log-likelihood of the im- 
age given <f> and a (assuming a gaussian noise model), 
while the second term may be interpreted as specifying 
a particular log-prior on a. That is, 
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with A = 2u 2 N j3. Thus, we may interpret E as being 
proportional to — log P(I, a\<f>), since 



P(J,a\<t>) = P(I\a,<t>)P{a) 



-i E (I,a\<t>) 
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(6) 
(7) 



How can we use this insight to improve our understand- 
ing of the algorithm? 

Under the maximum-likelihood approach, we would 
try to find the set of basis functions, <f>* , such that 



p(m 



arg max (log P(I\<f>)) 

4> 



P{I\a,<f>)P{a)da 



(8) 
(9) 



In other words, we are trying to find a set of <f>i that 
maximize the log-likelihood that the set of images could 
have arisen from a random process in which the <f>i are 
linearly mixed with statistically independent amplitudes 
distributed according to -^-e~^ s ^ ai > , with additive gaus- 
sian image noise. This is formally equivalent to mini- 
mizing the Kullback-Leibler (KL) distance between the 
actual joint probability of the images, P*(T), and our 



P(a) 



P(Ikt,4,) 



P(l,a\§) = P(l\a,§) P(a) 





Figure 1: Two-dimensional iso-probability plots of a, 
Cauchy prior, b, Gaussian likelihood, and c, their prod- 
uct. The axes on each plot are a\, 0,2- 

model of the joint probability based on independent 
causes, P(I\(f)), since 



KL[P*(I),P(I\<f>)] 



P*(I) 



= -Hp.-{logP(I\<j>)} (11) 

and Hp* = — f P* log P* is fixed, so maximizing 
(log P(I\(f))} minimizes KL. 

Unfortunately, all of this is easier said than done be- 
cause we have to integrate over the entire set of a; in 
Equation 9, which is computationally intractable. A 
reasonable approximation may be to assume that <7jv 
is small, in which case the dominant contribution to the 
integral is at the maximum of P(I, a\<f>). Thus, 

</>* = argmax (log[max_P(7|a, </>)P(a)] ) . (12) 

This is equivalent to the algorithm of Olshausen and 
Field (1996), as can be seen by comparing to Equation 3 
and using the definitions of Equations 4 and 5. The in- 
tuition for why this approximation works in practice is 
shown in Figure 1. The prior, P(a), is a product of 1-D 
"sparse" distributions, such as , } , , which are unimodal 

and peaked at zero. The likelihood, P(I\a, <f>), is a mul- 
tivariate gaussian, and since we are usually working in 
the overcomplete case (the number of basis functions ex- 
ceeds the dimensionality of the input) this will take the 
form of a gaussian ridge (or sandwich) that has its max- 
imum along the line (or plane, etc.) given by I = a<f>. 
The product of these two functions, P(I\a, <f)P[a), will 
have its maximum displaced away from the maximum 
along the gaussian ridge (i.e., away from the "perfect 
solution") and towards the origin, but also towards the 
ridges of the prior. Thus, the gradient with respect to <f> 
will tend to steer the gaussian ridge towards the ridges 
of the prior, which will in turn increase the volume of 
their product, or P(I\<f>). The reason we can get by with 
this approximation in this case is because we are working 
with a product of two fairly smooth, unimodal functions. 
If the functions were not so well behaved, then one can 
see that such an approximation might produce problems. 

4 Relation to Bell and Sejnowski 

Bell and Sejnowski (1995) describe an algorithm for "in- 
dependent component analysis" based on maximizing 



the mutual information between the inputs and out- 
puts of a neural network. Here, we show that this algo- 
rithm may be understood as solving the same maximum- 
likelihood problem described above (Section 3), except 
by making a different simplifying assumption. This has 
also been shown recently by Pearlmutter & Parra (1996) 
and Mackay (1996). 

Bell and Sejnowski examine the case where the num- 
ber of basis functions is equal to the number of inputs, 
and where the <f>i are linearly independent. In this case, 
there is a unique set of a; for which \I — a<j>\ 2 equals zero 
for any given image, I. In terms of the previous dis- 
cussion, P(I\a, <f>) is now a gaussian hump with a single 
maximum at a = I<f)~ l , rather than a gaussian ridge as 
in Figure lb. If we let <7jv go to zero in Equation 4, then 
P(I\a, <f>) becomes like a delta function and the integral 
of Equation 9 becomes 



p(m 
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P(I<t>~ v ) x I det ^ _1 | (14) 
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(16) 
By making the following definitions according to the con- 
vention of Bell and Sejnowski (1995), 

W = c/,- 1 (17) 

ui = Yfi-I (18) 

then, the gradient descent learning rule for W becomes 

cof ^ (19) 

detW • l yj 

This is precisely Bell and Sejnowski 's learning rule when 
the output non-linearity of their network, g[x), is equal 
to the cdf (cumulative density function) of the prior on 
the a,-, i.e., 



AWij oc -XS'(ui)Ij + 
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(20) 
(21) 



Thus, the independent component analysis algorithm 
of Bell and Sejnowski (1995) is formally equivalent to 
maximum likelihood in the case of no noise and a square 
system (dimensionality of output = dimensionality of in- 
put). It is easy to generalize this to the case when the 
number of outputs is less than the number of inputs, but 
not the other way around. When the number of outputs 
is greater than the effective dimensionality of the input 
($: of non-zero eigenvalues of the input covariance ma- 
trix), then the extra dimensions of the output will simply 
drop out. While this does not pose a problem for blind 
separation problems where the number of independent 
sources (dimensionality of a) is less than or equal to the 
number of mixed signals (dimensionality of I) , it will be- 
come a concern in the representation of images, where 
overcompleteness is a desirable feature (Simoncelli et al., 
1992). 



5 Lower-bound maximization 

A central idea behind the Helmholtz machine of Dayan 
et al. (1995), as well as the "mean field" theory of Saul 
et al. (1996), is the construction of an alternative proba- 
bility distribution, Q(a\T), that is used to obtain a lower- 
bound on logP(I\(f)). First, we rewrite logP(I\<f>) as 



logP(/|^) = log / Q(a\I) 



P( I\aA)P(a) 
Q(a\I) ' 



(22) 



Then, as long as Q is a probability (i.e., J Q = 1, Q > 0), 
we obtain by Jensen's inequality 



logP(J^) > 
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where Hq = — J Q(a\I) log Q(a\I)da. Thus, if we 
can construct Q(a\I) so that the integral is tractable, 
then we can do gradient ascent on a lower bound of 
(log P(I\(f))}. How good the bound is, though, depends 
on the Kullback-Leibler distance between Q(a\I) and 
P(a\I), or in other words, on how closely we can approx- 
imate P(a\I) with our tractable choice of Q. Typically, 
Q is chosen to be factorial, Q(a\I) = ILiqi(a,i\I), in which 
case 
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(26) 



where Hi = J q^a^a-ida-i and of = f qi(a,i)(a,i - m)' 2 da,i. 
Comparing Equation 26 to Equation 2, one can see 
that the sparse coding learning algorithm of Olshausen 
and Field (1996) effectively uses qi(ai) = 6(a; — Hi), with 
Hi chosen so as to minimize E (and hence maximize the 
lower bound of Equation 24). This choice would seem 
suboptimal, though, because we are getting zero entropy 
out of Hq (actually Hq = — oo, but we are ignoring the 
infinities here because it is the derivatives we really care 
about). If we could find a qi with higher entropy which 
also lowers the energy, then we could move the bound 
closer to the true log-likelihood. However, broadening 
qi (for example, by making it gaussian with adjustable 
Hi and <7 8 ) only affects the solution for h insofar as it 
low-pass filters the cost function, S, which has a simi- 
lar effect to simply lowering A. So, it is difficult to see 
that adding this extra complexity will improve matters. 
One apparent benefit of having non-zero o - ; is that there 
is now a growth-limiting term on the <f>i (second term 
on the right side of Eq. 26). Without such a term, the 
4>i will grow without bound, and so it is necessary in 
the algorithm of Olshausen and Field (1996) to keep the 
<f>i normalized (which is rather ad hoc by comparison). 
Preliminary investigation using a Gaussian qi and mini- 
mizing E with respect to both Hi an d v% for each image 
(but still keeping the <f>i normalized) does not reveal sig- 
nificant differences in the solution, but it deserves further 



study. It may also be worthwhile to try using a Q(a\I) 
that is defined by pairwise statistics (i.e., a covariance 
matrix on the a;). 

It should be noted that what is important here is 
the location of the maximum of whatever approximating 
function we use, not the absolute value of the bound per 
se. If the maximum of the lower-bound occurs at a sig- 
nificantly different point than the maximum of the true 
log-likelihood, then the approximation is not much help 
to us. 

6 Discussion 

What I think has been gained from this process is a bet- 
ter understanding of both the sparse coding algorithm of 
Olshausen and Field (1996) and the independent com- 
ponent analysis algorithm of Bell and Sejnowski (1995). 
Although neither of these algorithms was originally cast 
in maximum-likelihood terms, they are both essentially 
solving the same problem. The main difference between 
them is in the simplifying assumptions they make in 
order to deal with the intractable integration problem 
posed by Equation 9: Olshausen and Field's algorithm 
assumes low-noise (small cjv") and thus a peaky, uni- 
modal distribution on P(I, a\<f>) in order to justify eval- 
uating it at the maximum, whereas Bell and Sejnowski 
limit the dimensionality of the a; to equal the dimension- 
ality of the input and also assume no noise so that the 
integral becomes tractable. The maximum-likelihood 
framework also makes possible the link to techniques 
used in the Helmholtz machine (Dayan et al., 1995), 
which reveals that a better choice of approximating dis- 
tribution, Q, could potentially lead to improvements. 

A practical advantage of looking at the problem 
within this framework is that it suggests we could adapt 
the shape of the prior. For example, the prior on the 
a j need not be i.i.d., but could be shaped differently for 
each d{, e.g., P(cii) = ^ — e~' 3, ' s ( a, \ in order to best fit 

the data. Adapting /?; would be accomplished by letting 
it evolve along the gradient of (logP(I\(f))}. Using the 
approximation of Equation 12, this yields the learning 
rule: 

ft x -(S(ai))-—-^. (27) 

^ Pi opi 

A problem that may arise here, due to the fact that 
the full integral in Equation 9 is not being computed, is 
that there may be a bias toward non-informative flat pri- 
ors (since these will yield perfect reconstruction on each 
trial). An advantage of Bell and Sejnowski's algorithm 
in this case is that it essentially computes the full inte- 
gral in Equation 9 and so does not have this problem. 
For their algorithm, the maximum-likelihood framework 
prescribes a method for adapting the "generalized sig- 
moid" parameters p and r for shaping the prior (see pp. 
1137-8 of their paper), again by doing gradient ascent 
on the average log-likelihood. (See also Mackay, 1996, 
for other methods of parameterizing and adapting a fac- 
torial prior.) In cases where a statistically independent 
linear code may not be achieved (e.g., natural images), it 
may be advantageous to alter the prior so that informa- 
tion about pairwise or higher-order statistical dependen- 



cies among the a; may by incorporated into our model 
of P(a), for example using a Markov random field type 
model. 
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